/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. You can redistribute 
 * it and/or modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation, either version 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#include "mesh_io.h"

gctl::mesh_element::mesh_element()
{
    enabled = true;
    id = 0;
    type_tag = DEFAULT_INVALID_TAG;
}

gctl::mesh_element_group::mesh_element_group()
{
    enabled = true;
    type_tag = phys_group = geom_group = part_group = DEFAULT_INVALID_TAG;
}

gctl::mesh_data::mesh_data()
{
    enabled = true;
    d_type = NodeData;
}

void gctl::mesh_data::clear()
{
    enabled = true;
    str_tag.clear();
    real_tag.clear();
    int_tag.clear();
    vert_ptrs.clear();
    elem_ptrs.clear();
    val.clear();
    return;
}

gctl::mesh_physical_group::mesh_physical_group()
{
    enabled = true;
}


gctl::mesh_io::mesh_io()
{
    valid_node_size = 0;
    valid_elem_size = 0;
    valid_phys_size = 0;
    initialized_ = false;

    mesh_type_[0] = NoSet;
    mesh_type_[1] = _2NodeLine;
    mesh_type_[2] = _3NodeTriangle;
    mesh_type_[3] = _4NodeQuadrangle;
    mesh_type_[4] = _4NodeTetrahedron;
    mesh_type_[5] = _8NodeHexahedron;
    mesh_type_[6] = _6NodePrism;
    mesh_type_[7] = _5NodePyramid;
    mesh_type_[8] = _3NodeSecondOrderLine;
    mesh_type_[9] = _6NdoeSecondOrderLine;
    mesh_type_[10] = _9NodeSecondOrderQuadrangle;
    mesh_type_[11] = _10NodeSecondOrderTetrahedron;
    mesh_type_[12] = _27NodeSecondOrderHexahedron;
    mesh_type_[13] = _18NodeSecondOrderPrism;
    mesh_type_[14] = _14NodeSecondOrderPyramid;
    mesh_type_[15] = _1NodePoint;
    mesh_type_[16] = _8NodeSecondOrderQuadrangle;
    mesh_type_[17] = _20NdoeSecondOrderHexahedron;
    mesh_type_[18] = _15NodeSecondOrderPrism;
    mesh_type_[19] = _13NodeSecondOrderPyramid;
    mesh_type_[20] = _9NodeThirdOrderIncompleteTriangle;
    mesh_type_[21] = _10NdoeThirdOrderTriangle;
    mesh_type_[22] = _12NodeFourthOrderIncompleteTriangle;
    mesh_type_[23] = _15NodeFourthOrderTriangle;
    mesh_type_[24] = _15NodeFifthOrderCompleteTriangle;
    mesh_type_[25] = _21NodeFifthOrderCompleteTriangle;
    mesh_type_[26] = _4NodeThirdOrderEdge;
    mesh_type_[27] = _5NodeFourthOrderEdge;
    mesh_type_[28] = _6NodeFifthOrderEdge;
    mesh_type_[29] = _20NodeThirdOrderTetrahedron;
    mesh_type_[30] = _35NodeFourthOrderTetrahedron;
    mesh_type_[31] = _56NodeFifithOrderTetrahedron;
    mesh_type_[92] = _64NodeThirdOrderHexahedron;
    mesh_type_[93] = _125NodeFourthOrderHexahedron;

    mesh_type_name_[1] = "2-node line";
	mesh_type_name_[2] = "3-node triangle";
	mesh_type_name_[3] = "4-node quadrangle";
	mesh_type_name_[4] = "4-node tetrahedron";
	mesh_type_name_[5] = "8-node hexahedron";
	mesh_type_name_[6] = "6-node prism";
	mesh_type_name_[7] = "5-node pyramid";
	mesh_type_name_[8] = "3-node second order line";
	mesh_type_name_[9] = "6-ndoe second order line";
	mesh_type_name_[10] = "9-node second order quadrangle";
	mesh_type_name_[11] = "10-node second order tetrahedron";
	mesh_type_name_[12] = "27-node second order hexahedron";
	mesh_type_name_[13] = "18-node second order prism";
	mesh_type_name_[14] = "14-node second order pyramid";
	mesh_type_name_[15] = "1-node point";
	mesh_type_name_[16] = "8-node second order quadrangle";
	mesh_type_name_[17] = "20-ndoe second order hexahedron";
	mesh_type_name_[18] = "15-node second order prism";
	mesh_type_name_[19] = "13-node second order pyramid";
	mesh_type_name_[20] = "9-node third order incomplete triangle";
	mesh_type_name_[21] = "10-ndoe third order triangle";
	mesh_type_name_[22] = "12-node fourth order incomplete triangle";
	mesh_type_name_[23] = "15-node fourth order triangle";
	mesh_type_name_[24] = "15-node fifth order complete triangle";
	mesh_type_name_[25] = "21-node fifth order complete triangle";
	mesh_type_name_[26] = "4-node third order edge";
	mesh_type_name_[27] = "5-node fourth order edge";
	mesh_type_name_[28] = "6-node fifth order edge";
	mesh_type_name_[29] = "20-node third order tetrahedron";
	mesh_type_name_[30] = "35-node fourth order tetrahedron";
	mesh_type_name_[31] = "56-node fifith order tetrahedron";
	mesh_type_name_[92] = "64-node third order hexahedron";
	mesh_type_name_[93] = "125-node fourth order hexahedron";

	mesh_type_size_[1] = 2;
    mesh_type_size_[2] = 3;
    mesh_type_size_[3] = 4;
    mesh_type_size_[4] = 4;
    mesh_type_size_[5] = 8;
	mesh_type_size_[6] = 6;
    mesh_type_size_[7] = 5; 
    mesh_type_size_[8] = 3; 
    mesh_type_size_[9] = 6; 
    mesh_type_size_[10] = 6;
	mesh_type_size_[11] = 10; 
    mesh_type_size_[12] = 27; 
    mesh_type_size_[13] = 18; 
    mesh_type_size_[14] = 14; 
    mesh_type_size_[15] = 1;
	mesh_type_size_[16] = 8; 
    mesh_type_size_[17] = 20; 
    mesh_type_size_[18] = 15; 
    mesh_type_size_[19] = 13; 
    mesh_type_size_[20] = 9;
	mesh_type_size_[21] = 10; 
    mesh_type_size_[22] = 12; 
    mesh_type_size_[23] = 15; 
    mesh_type_size_[24] = 15; 
    mesh_type_size_[25] = 21;
	mesh_type_size_[26] = 4; 
    mesh_type_size_[27] = 5; 
    mesh_type_size_[28] = 6; 
    mesh_type_size_[29] = 20; 
    mesh_type_size_[30] = 35;
	mesh_type_size_[31] = 56; 
    mesh_type_size_[92] = 64; 
    mesh_type_size_[93] = 125;

    mesh2vtk_type_tag_[1] = 3;
    mesh2vtk_type_tag_[2] = 5;
    mesh2vtk_type_tag_[3] = 9;
    mesh2vtk_type_tag_[4] = 10;
    mesh2vtk_type_tag_[5] = 12;
    mesh2vtk_type_tag_[6] = 13;

    vtk2mesh_type_tag_[3] = 1;
    vtk2mesh_type_tag_[5] = 2;
    vtk2mesh_type_tag_[9] = 3;
    vtk2mesh_type_tag_[10] = 4;
    vtk2mesh_type_tag_[12] = 5;
    vtk2mesh_type_tag_[13] = 6;
}
    
gctl::mesh_io::~mesh_io()
{
    reset();
}

void gctl::mesh_io::reset()
{
    if (!nodes_.empty()) nodes_.clear();
    if (!elems_.empty()) elems_.clear();
    if (!phys_.empty()) phys_.clear();
    if (!nodes_tag_.empty()) nodes_tag_.clear();

    if (!groups_.empty())
    {
        for (size_t i = 0; i < groups_.size(); i++)
        {
            groups_[i].elem_ptrs.clear();
        }
        groups_.clear();
    }

    if (!datas_.empty())
    {
        for (size_t i = 0; i < datas_.size(); i++)
        {
            datas_[i].clear();
        }
        datas_.clear();
    }

    valid_node_size = 0;
    valid_elem_size = 0;
    valid_phys_size = 0;
    initialized_ = false;
}

void gctl::mesh_io::info(std::ostream &ss)
{
    int inttemp,temp_count;
	std::string strtemp;
	std::stringstream stemp;
	std::string iteam[4] = {"element type: ","physical group: ","geometrical group: ","mesh partition: "};

    ss << "number of nodes: " << valid_node_size << ", number of elements: " << valid_elem_size << std::endl;

    for (size_t g = 0; g < groups_.size(); g++)
    {
        if (!groups_[g].enabled) continue;

        ss << iteam[0] << mesh_type_name_[groups_[g].type_tag] << ", ";
        if (groups_[g].phys_group != DEFAULT_INVALID_TAG)
        {
            ss << iteam[1] << groups_[g].phys_group;
            if (!phys_.empty())
            {
                for (size_t p = 0; p < phys_.size(); p++)
                {
                    if (phys_[p].phys_tag == groups_[g].phys_group)
                    {
                        ss << " \"" << phys_[p].name << "\"";
                        break;
                    }
                }
            }
            ss << ", ";
        }
        if (groups_[g].geom_group != DEFAULT_INVALID_TAG) ss << iteam[2] << groups_[g].geom_group << ", ";
        if (groups_[g].part_group != DEFAULT_INVALID_TAG) ss << iteam[3] << groups_[g].part_group << ", ";
        ss << "number of elements: " << groups_[g].elem_ptrs.size() << std::endl;
    }

    size_t d_size;
    double min, max;
    for (size_t d = 0; d < datas_.size(); d++)
    {
        if (!datas_[d].enabled) continue;

        d_size = 0;
        min = GCTL_BDL_MAX;
        max = GCTL_BDL_MIN;

        if (datas_[d].d_type == NodeData)
        {
            ss << "nodedata: \"";
            for (size_t l = 0; l < datas_[d].val.size(); l++)
            {
                if (datas_[d].vert_ptrs[l]->id != DEFAULT_INVALID_TAG)
                {
                    min = std::min(min, datas_[d].val[l]);
                    max = std::max(max, datas_[d].val[l]);
                    d_size++;
                }
            }
        }
        else
        {
            ss << "elementdata: \"";
            for (size_t l = 0; l < datas_[d].val.size(); l++)
            {
                if (datas_[d].elem_ptrs[l]->enabled)
                {
                    min = std::min(min, datas_[d].val[l]);
                    max = std::max(max, datas_[d].val[l]);
                    d_size++;
                }
            }
        }

        if (d_size) ss << datas_[d].str_tag[0] << "\", data number: " << d_size << ", mini: " << min << ", maxi = " << max << std::endl;
        else ss << datas_[d].str_tag[0] << "\", no valid value found.\n";
    }
    
    return;
}

void gctl::mesh_io::enable_data(std::string dataname)
{
    if (dataname == "null")
    {
        for (size_t i = 0; i < datas_.size(); i++)
        {
            datas_[i].enabled = true;
        }
    }
    else
    {
        for (size_t i = 0; i < datas_.size(); i++)
        {
            if (datas_[i].str_tag.front() == dataname) datas_[i].enabled = true;
        }
    }
    return;
}

void gctl::mesh_io::disable_data(std::string dataname)
{
    if (dataname == "null")
    {
        for (size_t i = 0; i < datas_.size(); i++)
        {
            datas_[i].enabled = false;
        }
    }
    else
    {
        for (size_t i = 0; i < datas_.size(); i++)
        {
            if (datas_[i].str_tag.front() == dataname) datas_[i].enabled = false;
        }
    }
    return;
}

void gctl::mesh_io::enable_element(element_type_enum elem_type)
{
    if (elem_type == NoSet)
    {
        for (size_t e = 0; e < elems_.size(); e++)
        {
            elems_[e].enabled = true;
        }
    }
    else
    {
        for (size_t e = 0; e < elems_.size(); e++)
        {
            if (elems_[e].type == elem_type)
            {
                elems_[e].enabled = true;
            }
        }
    }

    update_indexing();
    return;
}

void gctl::mesh_io::enable_element(std::string phys_name)
{
    if (phys_.empty()) throw std::runtime_error("No physical groups found. From gctl::mesh_io::enable_element(...)");
    
    int p_tag;
    for (size_t p = 0; p < phys_.size(); p++)
    {
        if (phys_[p].name == phys_name)
        {
            p_tag = phys_[p].phys_tag;
        }
    }

    for (size_t e = 0; e < elems_.size(); e++)
    {
        if (elems_[e].int_tag[0] == p_tag) elems_[e].enabled = true;
    }

    update_indexing();
    return;
}

void gctl::mesh_io::enable_element(int tag, element_tag_enum tag_type)
{
    for (size_t e = 0; e < elems_.size(); e++)
    {
        if (tag_type == PhysicalTag && elems_[e].int_tag[0] == tag)
        {
            elems_[e].enabled = true;
        }

        if (tag_type == GeometryTag && elems_[e].int_tag[1] == tag)
        {
            elems_[e].enabled = true;
        }

        if (tag_type == PartitionTag && elems_[e].int_tag[2] == tag)
        {
            elems_[e].enabled = true;
        }
    }

    update_indexing();
    return;
}

void gctl::mesh_io::disable_element(element_type_enum elem_type)
{
    if (elem_type == NoSet)
    {
        for (size_t e = 0; e < elems_.size(); e++)
        {
            elems_[e].enabled = false;
        }
    }
    else
    {
        for (size_t e = 0; e < elems_.size(); e++)
        {
            if (elems_[e].type == elem_type)
            {
                elems_[e].enabled = false;
            }
        }
    }

    update_indexing();
    return;
}

void gctl::mesh_io::disable_element(std::string phys_name)
{
    if (phys_.empty()) throw std::runtime_error("No physical groups found. From gctl::mesh_io::disable_element(...)");
    
    int p_tag;
    for (size_t p = 0; p < phys_.size(); p++)
    {
        if (phys_[p].name == phys_name)
        {
            p_tag = phys_[p].phys_tag;
        }
    }

    for (size_t e = 0; e < elems_.size(); e++)
    {
        if (elems_[e].int_tag[0] == p_tag) elems_[e].enabled = false;
    }

    update_indexing();
    return;
}

void gctl::mesh_io::disable_element(int tag, element_tag_enum tag_type)
{
    for (size_t e = 0; e < elems_.size(); e++)
    {
        if (tag_type == PhysicalTag && elems_[e].int_tag[0] == tag)
        {
            elems_[e].enabled = false;
        }

        if (tag_type == GeometryTag && elems_[e].int_tag[1] == tag)
        {
            elems_[e].enabled = false;
        }

        if (tag_type == PartitionTag && elems_[e].int_tag[2] == tag)
        {
            elems_[e].enabled = false;
        }
    }

    update_indexing();
    return;
}

const gctl::array<int> &gctl::mesh_io::get_node_tag()
{
    return nodes_tag_;
}

const gctl::array<gctl::vertex3dc> &gctl::mesh_io::get_nodes()
{
    return nodes_;
}

const gctl::array<gctl::mesh_element> &gctl::mesh_io::get_elems()
{
    return elems_;
}

void gctl::mesh_io::export_tags_to_data(element_tag_enum tag_type)
{
    mesh_data tmp_data;
    if (tag_type == NodeTag && (!nodes_tag_.empty()))
    {
        tmp_data.enabled = true;
        tmp_data.d_type = NodeData;
        tmp_data.str_tag.resize(1, "Node Tag");
        tmp_data.real_tag.resize(1, 0.0);
        tmp_data.int_tag.resize(3, 0);
        tmp_data.int_tag[1] = 1;
        tmp_data.int_tag[2] = valid_node_size;

        tmp_data.val.resize(valid_node_size);
        tmp_data.vert_ptrs.resize(valid_node_size);

        size_t c = 0;
        for (size_t i = 0; i < nodes_.size(); i++)
        {
            if (nodes_[i].id != DEFAULT_INVALID_TAG)
            {
                tmp_data.val[c] = (double) nodes_tag_[i];
                tmp_data.vert_ptrs[c] = nodes_.get(i);
                c++;
            }
        }

        datas_.push_back(tmp_data);
    }
    else
    {
        int t = 0;
        tmp_data.enabled = true;
        tmp_data.d_type = ElemData;

        if (tag_type == PhysicalTag) {tmp_data.str_tag.resize(1, "Physical Tag"); t = 0;}
        else if (tag_type == GeometryTag) {tmp_data.str_tag.resize(1, "Geometry Tag"); t = 1;}
        else if (tag_type == PartitionTag) {tmp_data.str_tag.resize(1, "Partition Tag"); t = 2;}

        tmp_data.real_tag.resize(1, 0.0);
        tmp_data.int_tag.resize(3, 0);
        tmp_data.int_tag[1] = 1;
        tmp_data.int_tag[2] = valid_elem_size;

        tmp_data.val.resize(valid_elem_size);
        tmp_data.elem_ptrs.resize(valid_elem_size);

        size_t c = 0;
        for (size_t i = 0; i < elems_.size(); i++)
        {
            if (elems_[i].enabled)
            {
                tmp_data.val[c] = (double) elems_[i].int_tag[t];
                tmp_data.elem_ptrs[c] = elems_.get(i);
                c++;
            }
        }

        datas_.push_back(tmp_data);
    }
    return;
}

void gctl::mesh_io::export_elements_to(array<triangle> &tris)
{
    size_t s = 0;
    for (size_t i = 0; i < elems_.size(); i++)
    {
        if (elems_[i].type_tag == 2) s++;
    }
    
    tris.resize(s);

    s = 0;
    for (size_t i = 0; i < elems_.size(); i++)
    {
        if (elems_[i].type_tag == 2)
        {
            tris[s].id = elems_[i].id;
            tris[s].vert[0] = elems_[i].vert_ptrs[0];
            tris[s].vert[1] = elems_[i].vert_ptrs[1];
            tris[s].vert[2] = elems_[i].vert_ptrs[2];
            s++;
        }
    }
    return;
}

void gctl::mesh_io::export_elements_to(array<tetrahedron> &tets)
{
    size_t s = 0;
    for (size_t i = 0; i < elems_.size(); i++)
    {
        if (elems_[i].type_tag == 4) s++;
    }
    
    tets.resize(s);

    s = 0;
    for (size_t i = 0; i < elems_.size(); i++)
    {
        if (elems_[i].type_tag == 4)
        {
            tets[s].set(elems_[i].vert_ptrs[0], elems_[i].vert_ptrs[1], 
                elems_[i].vert_ptrs[2], elems_[i].vert_ptrs[3], elems_[i].id);
            s++;
        }
    }
    return;
}

void gctl::mesh_io::read_triangle_ascii(std::string filename, index_packed_e is_packed)
{
    array<vertex2dc> node_2d;
    array<int> node_tag;
    read_Triangle_node(filename, node_2d, is_packed, &node_tag);
    if (!node_tag.empty()) nodes_tag_ = node_tag;

    array<triangle2d> tri_2d;
    matrix<double> tri_attri;
    read_Triangle_element(filename, tri_2d, node_2d, is_packed, &tri_attri);

    valid_node_size = node_2d.size();
    valid_elem_size = tri_2d.size();
    valid_phys_size = 0;

    nodes_.resize(valid_node_size);
    for (size_t i = 0; i < valid_node_size; i++)
    {
        nodes_[i].id = i;
        nodes_[i].x = node_2d[i].x;
        nodes_[i].y = node_2d[i].y;
        nodes_[i].z = 0.0;
    }

    elems_.resize(valid_elem_size);
    for (size_t i = 0; i < valid_elem_size; i++)
    {
        elems_[i].enabled = true;
        elems_[i].id = i;
        elems_[i].int_tag.resize(3, DEFAULT_INVALID_TAG);
        elems_[i].neigh_ptrs.resize(3, nullptr);
        elems_[i].type = _3NodeTriangle;
        elems_[i].type_tag = 2;

        elems_[i].vert_ptrs.resize(3);
        elems_[i].vert_ptrs[0] = nodes_.get(tri_2d[i].vert[0]->id);
        elems_[i].vert_ptrs[1] = nodes_.get(tri_2d[i].vert[1]->id);
        elems_[i].vert_ptrs[2] = nodes_.get(tri_2d[i].vert[2]->id);
    }

    if (!tri_attri.empty())
    {
        for (size_t i = 0; i < valid_elem_size; i++)
        {
            elems_[i].int_tag[0] = 0;
            elems_[i].int_tag[1] = (int) tri_attri[i][0];
        }
    }
    
    // 整理单元体组
    mesh_element_group tmp_group;
    tmp_group.type_tag = elems_[0].type_tag;
    tmp_group.phys_group = elems_[0].int_tag[0];
    tmp_group.geom_group = elems_[0].int_tag[1];
    tmp_group.part_group = elems_[0].int_tag[2];
    tmp_group.elem_ptrs.push_back(elems_.get(0));
    groups_.push_back(tmp_group);

    bool not_found;
    for (size_t i = 1; i < elems_.size(); i++)
    {
        not_found = true;
        for (size_t g = 0; g < groups_.size(); g++)
        {
            if (groups_[g].type_tag == elems_[i].type_tag &&
                groups_[g].phys_group == elems_[i].int_tag[0] &&
                groups_[g].geom_group == elems_[i].int_tag[1] &&
                groups_[g].part_group == elems_[i].int_tag[2])
            {
                groups_[g].elem_ptrs.push_back(elems_.get(i));
                not_found = false;
                break;
            }
        }
        
        if (not_found)
        {
            tmp_group.elem_ptrs.clear();
            tmp_group.type_tag = elems_[i].type_tag;
            tmp_group.phys_group = elems_[i].int_tag[0];
            tmp_group.geom_group = elems_[i].int_tag[1];
            tmp_group.part_group = elems_[i].int_tag[2];
            tmp_group.elem_ptrs.push_back(elems_.get(i));
            groups_.push_back(tmp_group);
        }
    }

    std::string n_file = filename + ".neigh";
    if (!access(n_file.c_str(), R_OK))
    {
        read_Triangle_neighbor(filename, tri_2d, is_packed);

        for (size_t i = 0; i < valid_elem_size; i++)
        {
            for (size_t j = 0; j < 3; j++)
            {
                if (tri_2d[i].neigh[j] != nullptr) elems_[i].neigh_ptrs[j] = elems_.get(tri_2d[i].neigh[j]->id);
            }
        }
    }

    initialized_ = true;
    return;
}

void gctl::mesh_io::read_gmsh_v2_ascii(std::string filename, index_packed_e is_packed)
{
    if (initialized_ == true)
    {
        throw std::runtime_error("The object is already initialized. From gctl::mesh_io::read_mesh(...)");
    }

    std::ifstream infile;
    open_infile(infile, filename, ".msh");

    std::string tmp_str, tmp_str2;
    std::stringstream tmp_ss;

    // 检查mesh文件版本
    while(getline(infile,tmp_str))
    {
        if (tmp_str == "$MeshFormat")
        {
            getline(infile, tmp_str);
            if (tmp_str != "2.2 0 8")
            {
                throw std::runtime_error("Unsupported mesh file format. From gctl::mesh_io::read_mesh(...)");
            }
            break;
        }
    }

    //读入模型空间顶点集和物理分组
    size_t p_size = 0;
    size_t n_size = 0;
    
    while(getline(infile,tmp_str))
    {
        if (tmp_str == "$PhysicalNames")
        {
            getline(infile, tmp_str);
            gctl::str2ss(tmp_str, tmp_ss);
            tmp_ss >> p_size;
            valid_phys_size = p_size;

            phys_.resize(p_size);
            for (size_t i = 0; i < p_size; i++)
            {
                getline(infile, tmp_str);
                str2ss(tmp_str, tmp_ss);
                tmp_ss >> phys_[i].dim_tag >> phys_[i].phys_tag >> tmp_str2;
                // 去掉两端的双引号
                replace_all(phys_[i].name, tmp_str2, "\"", "");
            }
        }

        if (tmp_str == "$Nodes")
        {
            getline(infile, tmp_str);
            str2ss(tmp_str, tmp_ss);
            tmp_ss >> n_size; //第一个数为顶点个数
            nodes_.resize(n_size);
            valid_node_size = n_size;

            for (int i = 0; i < n_size; i++)
            {
                getline(infile, tmp_str);
                str2ss(tmp_str, tmp_ss);
                tmp_ss >> nodes_[i].id >> nodes_[i].x >> nodes_[i].y >> nodes_[i].z;
                if (is_packed == NotPacked) nodes_[i].id -= 1;
            }
            break;
        }
    }

    // 读入模型空间元素集
    int i_size, attri_num, vt_idx;
    while(getline(infile,tmp_str))
    {
        if (tmp_str == "$Elements")
        {
            getline(infile,tmp_str);
            gctl::str2ss(tmp_str, tmp_ss);
            tmp_ss >> i_size; //第一个数为元素个数
            elems_.resize(i_size);
            valid_elem_size = i_size;

            for (size_t i = 0; i < i_size; i++)
            {
                getline(infile,tmp_str);
                str2ss(tmp_str, tmp_ss);
                tmp_ss >> elems_[i].id >> elems_[i].type_tag >> attri_num;

                elems_[i].type = mesh_type_[elems_[i].type_tag];
                elems_[i].vert_ptrs.resize(mesh_type_size_[elems_[i].type_tag]);

                if (attri_num >= 3) // default tags will be assgined to DEFAULT_INVALID_TAG
                {
                    elems_[i].int_tag.resize(attri_num);
                    for (size_t a = 0; a < attri_num; a++)
                    {
                        tmp_ss >> elems_[i].int_tag[a];
                    }
                }
                else if (attri_num == 2)
                {
                    elems_[i].int_tag.resize(3, DEFAULT_INVALID_TAG);
                    tmp_ss >> elems_[i].int_tag[0];
                    tmp_ss >> elems_[i].int_tag[1];
                }
                else if (attri_num == 1)
                {
                    elems_[i].int_tag.resize(3, DEFAULT_INVALID_TAG);
                    tmp_ss >> elems_[i].int_tag[0];
                }
                else
                {
                    elems_[i].int_tag.resize(3, DEFAULT_INVALID_TAG);
                }

                for (size_t v = 0; v < elems_[i].vert_ptrs.size(); v++)
                {
                    tmp_ss >> vt_idx;
                    if (is_packed == Packed) elems_[i].vert_ptrs[v] = nodes_.get(vt_idx);
                    else elems_[i].vert_ptrs[v] = nodes_.get(vt_idx - 1);
                }
            }
            break;
        }
    }

    // 读入数据模块
    mesh_data tmp_data;
    while(getline(infile,tmp_str))
    {
        if (!infile.good()) break;

        if (tmp_str == "$NodeData" || tmp_str == "$ElementData")
        {
            tmp_data.clear();

            if (tmp_str == "$NodeData") tmp_data.d_type = NodeData;
            if (tmp_str == "$ElementData") tmp_data.d_type = ElemData;

            getline(infile,tmp_str);
            str2ss(tmp_str, tmp_ss);
            tmp_ss >> i_size;
            tmp_data.str_tag.resize(i_size);
            for (size_t i = 0; i < i_size; i++)
            {
                getline(infile, tmp_str);
                replace_all(tmp_str2, tmp_str, "\"", "");
                tmp_data.str_tag[i] = tmp_str2;
            }

            getline(infile,tmp_str);
            str2ss(tmp_str, tmp_ss);
            tmp_ss >> i_size;
            tmp_data.real_tag.resize(i_size);
            for (size_t i = 0; i < i_size; i++)
            {
                getline(infile,tmp_str);
                str2ss(tmp_str, tmp_ss);
                tmp_ss >> tmp_data.real_tag[i];
            }

            getline(infile,tmp_str);
            str2ss(tmp_str, tmp_ss);
            tmp_ss >> i_size;
            tmp_data.int_tag.resize(i_size);
            for (size_t i = 0; i < i_size; i++)
            {
                getline(infile,tmp_str);
                str2ss(tmp_str, tmp_ss);
                tmp_ss >> tmp_data.int_tag[i];
            }

            if (tmp_data.d_type == NodeData)
            {
                tmp_data.vert_ptrs.resize(tmp_data.int_tag.back());
                tmp_data.val.resize(tmp_data.int_tag.back());
                for (size_t i = 0; i < tmp_data.val.size(); i++)
                {
                    getline(infile,tmp_str);
                    str2ss(tmp_str, tmp_ss);
                    tmp_ss >> vt_idx >> tmp_data.val[i];
                    if (is_packed == Packed) tmp_data.vert_ptrs[i] = nodes_.get(vt_idx);
                    else tmp_data.vert_ptrs[i] = nodes_.get(vt_idx - 1);
                }
            }
            else
            {
                tmp_data.elem_ptrs.resize(tmp_data.int_tag.back());
                tmp_data.val.resize(tmp_data.int_tag.back());
                for (size_t i = 0; i < tmp_data.val.size(); i++)
                {
                    getline(infile,tmp_str);
                    str2ss(tmp_str, tmp_ss);
                    tmp_ss >> vt_idx >> tmp_data.val[i];
                    if (is_packed == Packed) tmp_data.elem_ptrs[i] = elems_.get(vt_idx);
                    else tmp_data.elem_ptrs[i] = elems_.get(vt_idx - 1);
                }
            }

            datas_.push_back(tmp_data);
            continue;
        }
    }
    infile.close();

    // 整理单元体组
    mesh_element_group tmp_group;
    tmp_group.type_tag = elems_[0].type_tag;
    tmp_group.phys_group = elems_[0].int_tag[0];
    tmp_group.geom_group = elems_[0].int_tag[1];
    tmp_group.part_group = elems_[0].int_tag[2];
    tmp_group.elem_ptrs.push_back(elems_.get(0));
    groups_.push_back(tmp_group);

    bool not_found;
    for (size_t i = 1; i < elems_.size(); i++)
    {
        not_found = true;
        for (size_t g = 0; g < groups_.size(); g++)
        {
            if (groups_[g].type_tag == elems_[i].type_tag &&
                groups_[g].phys_group == elems_[i].int_tag[0] &&
                groups_[g].geom_group == elems_[i].int_tag[1] &&
                groups_[g].part_group == elems_[i].int_tag[2])
            {
                groups_[g].elem_ptrs.push_back(elems_.get(i));
                not_found = false;
                break;
            }
        }
        
        if (not_found)
        {
            tmp_group.elem_ptrs.clear();
            tmp_group.type_tag = elems_[i].type_tag;
            tmp_group.phys_group = elems_[i].int_tag[0];
            tmp_group.geom_group = elems_[i].int_tag[1];
            tmp_group.part_group = elems_[i].int_tag[2];
            tmp_group.elem_ptrs.push_back(elems_.get(i));
            groups_.push_back(tmp_group);
        }
    }

    initialized_ = true;
    return;
}

void gctl::mesh_io::save_gmsh_v2_ascii(std::string filename, index_packed_e is_packed)
{
    if (initialized_ == false)
    {
        throw std::runtime_error("The object is not initialized. From gctl::mesh_io::save_mesh(...)");
    }

    std::ofstream outfile;
    open_outfile(outfile, filename, ".msh");

    outfile << "$MeshFormat\n2.2 0 8\n$EndMeshFormat\n";

    if (valid_phys_size)
    {
        outfile << "$PhysicalNames\n" << valid_phys_size << "\n";
        for (size_t i = 0; i < phys_.size(); i++)
        {
            if (phys_[i].enabled) outfile << phys_[i].dim_tag << " " << phys_[i].phys_tag << " \"" << phys_[i].name << "\"\n";
        }
        outfile << "$EndPhysicalNames\n";
    }

    outfile << "$Nodes\n" << valid_node_size << std::endl;
    for (size_t i = 0; i < nodes_.size(); i++)
    {
        if (nodes_[i].id == DEFAULT_INVALID_TAG) continue;
    
        if (is_packed == NotPacked) outfile << nodes_[i].id + 1 << " ";
        else outfile << nodes_[i].id << " ";

        outfile << std::setprecision(16) << nodes_[i].x << " " << nodes_[i].y << " " << nodes_[i].z << std::endl;
    }
    outfile<<"$EndNodes\n";

    size_t a_size;
    outfile << "$Elements\n" << valid_elem_size << std::endl;
    for (size_t i = 0; i < elems_.size(); i++)
    {
        if (!elems_[i].enabled) continue;

        if (is_packed == NotPacked) outfile << elems_[i].id+1 << " " << elems_[i].type_tag << " ";
        else outfile << elems_[i].id << " " << elems_[i].type_tag << " ";

        a_size = 0;
        for (size_t a = 0; a < elems_[i].int_tag.size(); a++)
        {
            if (elems_[i].int_tag[a] != DEFAULT_INVALID_TAG) a_size++;
        }

        outfile << a_size;
        for (size_t a = 0; a < elems_[i].int_tag.size(); a++)
        {
            if (elems_[i].int_tag[a] != DEFAULT_INVALID_TAG) outfile << " " << elems_[i].int_tag[a];
        }

        for (size_t v = 0; v < elems_[i].vert_ptrs.size(); v++)
        {
            if (is_packed == NotPacked) outfile << " " << elems_[i].vert_ptrs[v]->id + 1;
            else outfile << " " << elems_[i].vert_ptrs[v]->id;
        }
        outfile << std::endl;
    }
    outfile << "$EndElements\n";

    if (!datas_.empty())
    {
        size_t d_size;
        for (size_t i = 0; i < datas_.size(); i++)
        {
            if (!datas_[i].enabled) continue;

            if (datas_[i].d_type == NodeData)
            {
                d_size = 0;
                for (size_t n = 0; n < datas_[i].val.size(); n++)
                {
                    if (datas_[i].vert_ptrs[n]->id != DEFAULT_INVALID_TAG) d_size++;
                }
                
                if (d_size)
                {
                    outfile << "$NodeData\n";
                    outfile << datas_[i].str_tag.size() << std::endl;
                    for (size_t a = 0; a < datas_[i].str_tag.size(); a++)
                    {
                        outfile << "\"" << datas_[i].str_tag[a] << "\"\n";
                    }

                    outfile << datas_[i].real_tag.size() << std::endl;
                    for (size_t a = 0; a < datas_[i].real_tag.size(); a++)
                    {
                        outfile << datas_[i].real_tag[a] << "\n";
                    }

                    outfile << datas_[i].int_tag.size() << std::endl;
                    for (size_t a = 0; a < datas_[i].int_tag.size() - 1; a++)
                    {
                        outfile << datas_[i].int_tag[a] << "\n";
                    }
                    outfile << d_size << "\n";

                    for (size_t a = 0; a < datas_[i].val.size(); a++)
                    {
                        if (datas_[i].vert_ptrs[a]->id != DEFAULT_INVALID_TAG && is_packed == Packed) outfile << datas_[i].vert_ptrs[a]->id << " " << std::setprecision(12) <<  datas_[i].val[a] << "\n";
                        else if (datas_[i].vert_ptrs[a]->id != DEFAULT_INVALID_TAG) outfile << datas_[i].vert_ptrs[a]->id + 1 << " " << std::setprecision(12) <<  datas_[i].val[a] << "\n";
                    }
                    outfile << "$EndNodeData\n";
                }
            }
            else
            {
                d_size = 0;
                for (size_t n = 0; n < datas_[i].val.size(); n++)
                {
                    if (datas_[i].elem_ptrs[n]->enabled) d_size++;
                }
                
                if (d_size)
                {
                    outfile << "$ElementData\n";
                    outfile << datas_[i].str_tag.size() << std::endl;
                    for (size_t a = 0; a < datas_[i].str_tag.size(); a++)
                    {
                        outfile << "\"" << datas_[i].str_tag[a] << "\"\n";
                    }

                    outfile << datas_[i].real_tag.size() << std::endl;
                    for (size_t a = 0; a < datas_[i].real_tag.size(); a++)
                    {
                        outfile << datas_[i].real_tag[a] << "\n";
                    }

                    outfile << datas_[i].int_tag.size() << std::endl;
                    for (size_t a = 0; a < datas_[i].int_tag.size() - 1; a++)
                    {
                        outfile << datas_[i].int_tag[a] << "\n";
                    }
                    outfile << d_size << "\n";

                    for (size_t a = 0; a < datas_[i].val.size(); a++)
                    {
                        if (datas_[i].elem_ptrs[a]->enabled && is_packed == Packed) outfile << datas_[i].elem_ptrs[a]->id << " " << std::setprecision(12) <<  datas_[i].val[a] << "\n";
                        else if (datas_[i].elem_ptrs[a]->enabled) outfile << datas_[i].elem_ptrs[a]->id + 1 << " " << std::setprecision(12) <<  datas_[i].val[a] << "\n";
                    }
                    outfile << "$EndElementData\n";
                }
            }
        }
    }

    outfile.close();
    return;
}

void gctl::mesh_io::save_vtk_legacy_ascii(std::string filename)
{
    if (initialized_ == false)
    {
        throw std::runtime_error("The object is not initialized. From gctl::mesh_io::export_vtk(...)");
    }

    std::ofstream outfile;
    open_outfile(outfile, filename, ".vtk");

    outfile << "# vtk DataFile Version 2.0\nGenerated by gctl::mesh_io\nASCII\nDATASET UNSTRUCTURED_GRID\nPOINTS " << valid_node_size <<" FLOAT\n";
	for (int i = 0; i < nodes_.size(); i++)
	{
        if (nodes_[i].id == DEFAULT_INVALID_TAG) continue;

		outfile << std::setprecision(16) << nodes_[i].x << " " << nodes_[i].y << " " << nodes_[i].z << std::endl;
	}

	//计算一下CELLS的总长 注意类型名算一个 所以要加1
    size_t totalCellNum = 0;
	for (int i = 0; i < elems_.size(); i++)
	{
        if (elems_[i].enabled) totalCellNum += mesh_type_size_[elems_[i].type_tag] + 1;
	}
	outfile << "CELLS " << valid_elem_size << " " << totalCellNum << std::endl;

	for (int i = 0; i < elems_.size(); i++)
	{
        if (!elems_[i].enabled) continue;

		outfile << mesh_type_size_[elems_[i].type_tag] << " ";
		for (int j = 0; j < elems_[i].vert_ptrs.size(); j++)
		{
			outfile << elems_[i].vert_ptrs[j]->id << " ";
		}
		outfile << std::endl;
	}

	outfile << "CELL_TYPES " << valid_elem_size << std::endl;
	for (int i = 0; i < elems_.size(); i++)
	{
		if (elems_[i].enabled) outfile << mesh2vtk_type_tag_[elems_[i].type_tag] << "\n";
	}

    bool nd_head = true, ed_head = true;
    size_t d_size;
    std::string temp_dataName;
    if (!datas_.empty())
    {
        for (size_t i = 0; i < datas_.size(); i++)
        {
            if (!datas_[i].enabled) continue;

            if (datas_[i].d_type == NodeData)
            {
                d_size = 0;
                for (size_t n = 0; n < datas_[i].val.size(); n++)
                {
                    if (datas_[i].vert_ptrs[n]->id != DEFAULT_INVALID_TAG) d_size++;
                }
                
                if (d_size && d_size == valid_node_size)
                {
                    if (nd_head) {outfile<<"POINT_DATA "<< d_size << std::endl; nd_head = false; ed_head = true;}

                    gctl::replace_all(temp_dataName, datas_[i].str_tag.front()," ","_");
                    outfile << "SCALARS " << temp_dataName << " FLOAT\nLOOKUP_TABLE default\n";

                    for (size_t a = 0; a < datas_[i].val.size(); a++)
                    {
                        if (datas_[i].vert_ptrs[a]->id != DEFAULT_INVALID_TAG) outfile << std::setprecision(16) <<  datas_[i].val[a] << "\n";
                    }
                }
            }
            else
            {
                d_size = 0;
                for (size_t n = 0; n < datas_[i].val.size(); n++)
                {
                    if (datas_[i].elem_ptrs[n]->enabled) d_size++;
                }
                
                if (d_size && d_size == valid_elem_size)
                {
                    if (ed_head) {outfile<<"CELL_DATA "<< d_size << std::endl; ed_head = false; nd_head = true;}

                    gctl::replace_all(temp_dataName, datas_[i].str_tag.front()," ","_");
                    outfile << "SCALARS " << temp_dataName << " FLOAT\nLOOKUP_TABLE default\n";

                    for (size_t a = 0; a < datas_[i].val.size(); a++)
                    {
                        if (datas_[i].elem_ptrs[a]->enabled) outfile << std::setprecision(16) <<  datas_[i].val[a] << "\n";
                    }
                }
            }
        }
    }

    outfile.close();
    return;
}

void gctl::mesh_io::export_data_to_lld(std::string filename, std::string dataname, double refr, double refR)
{
    if (dataname == "null")
    {
        int tmp_size;
        std::string tmp_name;
        point3ds ps;
        point3dc pc;

        std::ofstream ofile;
        for (size_t i = 0; i < datas_.size(); i++)
        {
            gctl::replace_all(tmp_name, datas_[i].str_tag.front()," ","_");
            gctl::replace_all(tmp_name, tmp_name,"/","");
            open_outfile(ofile, filename + "_" + tmp_name, ".txt");

            if (datas_[i].d_type == NodeData)
            {
                for (size_t n = 0; n < datas_[i].vert_ptrs.size(); n++)
                {
                    if (datas_[i].vert_ptrs[n]->id != DEFAULT_INVALID_TAG)
                    {
                        ps = datas_[i].vert_ptrs[n]->c2s();
                        ps.rad -= ellipse_radius_2d(refR, refr, ps.lat*M_PI/180.0);

                        ofile << ps.lon << " " << ps.lat << " " << ps.rad << " " << datas_[i].val[n] << "\n";
                    }
                }
                
            }
            else
            {
                for (size_t e = 0; e < valid_elem_size; e++)
                {
                    if (elems_[e].enabled)
                    {
                        pc = point3dc(0.0, 0.0, 0.0);
                        for (size_t v = 0; v < elems_[e].vert_ptrs.size(); v++)
                        {
                            pc.x += elems_[e].vert_ptrs[v]->x;
                            pc.y += elems_[e].vert_ptrs[v]->y;
                            pc.z += elems_[e].vert_ptrs[v]->z;
                        }
                        pc.x /= elems_[e].vert_ptrs.size();
                        pc.y /= elems_[e].vert_ptrs.size();
                        pc.z /= elems_[e].vert_ptrs.size();

                        ps = pc.c2s();
                        ps.rad -= ellipse_radius_2d(refR, refr, ps.lat*M_PI/180.0);

                        ofile << ps.lon << " " << ps.lat << " " << ps.rad << " " << datas_[i].val[e] << "\n";
                    }
                }

                for (size_t e = 0; e < datas_[i].elem_ptrs.size(); e++)
                {
                    if (datas_[i].elem_ptrs[e]->enabled)
                    {
                        pc = point3dc(0.0, 0.0, 0.0);
                        tmp_size = datas_[i].elem_ptrs[e]->vert_ptrs.size();
                        for (size_t v = 0; v < datas_[i].elem_ptrs[e]->vert_ptrs.size(); v++)
                        {
                            pc.x += datas_[i].elem_ptrs[e]->vert_ptrs[v]->x;
                            pc.y += datas_[i].elem_ptrs[e]->vert_ptrs[v]->y;
                            pc.z += datas_[i].elem_ptrs[e]->vert_ptrs[v]->z;
                        }
                        pc.x /= tmp_size;
                        pc.y /= tmp_size;
                        pc.z /= tmp_size;

                        ps = pc.c2s();
                        ps.rad -= ellipse_radius_2d(refR, refr, ps.lat*M_PI/180.0);

                        ofile << ps.lon << " " << ps.lat << " " << ps.rad << " " << datas_[i].val[e] << "\n";
                    }
                }
            }
            ofile.close();
        }
    }
    else
    {

    }
    return;
}

void gctl::mesh_io::update_indexing()
{
    for (size_t n = 0; n < nodes_.size(); n++)
    {
        nodes_[n].id = DEFAULT_INVALID_TAG;
    }
    
    int new_id = 0;
    for (size_t e = 0; e < elems_.size(); e++)
    {
        if (elems_[e].enabled)
        {
            elems_[e].id = new_id; new_id++;
            for (size_t n = 0; n < elems_[e].vert_ptrs.size(); n++)
            {
                elems_[e].vert_ptrs[n]->id = 1;
            }
        }
    }
    valid_elem_size = new_id;

    new_id = 0;
    for (size_t n = 0; n < nodes_.size(); n++)
    {
        if (nodes_[n].id == 1)
        {
            nodes_[n].id = new_id; new_id++;
        }
    }
    valid_node_size = new_id;

    valid_phys_size = 0;
    for (size_t p = 0; p < phys_.size(); p++)
    {
        phys_[p].enabled = false;

        for (size_t e = 0; e < elems_.size(); e++)
        {
            if (elems_[e].enabled && elems_[e].int_tag[0] == phys_[p].phys_tag)
            {
                phys_[p].enabled = true;
            }
        }

        if (phys_[p].enabled) valid_phys_size++;
    }

    for (size_t g = 0; g < groups_.size(); g++)
    {
        groups_[g].enabled = false;

        for (size_t e = 0; e < groups_[g].elem_ptrs.size(); e++)
        {
            if (groups_[g].elem_ptrs[e]->enabled)
            {
                groups_[g].enabled = true;
                break;
            }
        }
    }
    return;
}